Seurat4.0系列教程4：整合分析

分享是一种态度

scRNA-seq整合简介

对两个或两个以上单细胞数据集的整合分析提出了独特的挑战。特别是，在标准工作流下，识别存在于多个数据集中的基因可能存在问题。Seurat v4 包括一组方法，以匹配（或"对齐"）跨数据集共同的基因。这些方法首先识别处于匹配生物状态的交叉数据集对（"锚点"），既可用于纠正数据集之间的技术差异（即批次效应校正），又可用于对整个实验条件进行比较scRNA-seq分析。

下面，我们演示ScRNA-seq 整合的方法，在ctrl或干扰素刺激状态^[1]下对人体免疫细胞 （PBMC） 进行比较分析。

整合目标

以下教程旨在为您概述使用 Seurat 集成程序可能的复杂细胞类型的比较分析。在这里，我们讨论几个关键目标：

• 创建"整合"数据，用于下游分析
• 识别两个数据集中存在的细胞类型
• 获取在对照和刺激细胞中保守的细胞类型标记
• 比较数据集，找到细胞类型对刺激的特定反应

设置seurat对象

为了方便起见，我们通过我们的SeuratData^[2]包分发此数据集。

library(Seurat)
library(SeuratData)
library(patchwork)

# install dataset
InstallData("ifnb")

# load dataset
LoadData("ifnb")

# split the dataset into a list of two seurat objects (stim and CTRL)
ifnb.list <- SplitObject(ifnb, split.by = "stim")

# normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = ifnb.list)

执行整合

然后，我们使用该功能识别锚点，以 Seurat 对象作为输入，并使用这些锚点将两个数据集整合在一起。

immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features)

# this command creates an 'integrated' data assay
immune.combined <- IntegrateData(anchorset = immune.anchors)

执行整合分析

现在，我们可以对所有细胞进行单次整合分析！

# specify that we will perform downstream analysis on the corrected data note that the original
# unmodified data still resides in the 'RNA' assay
DefaultAssay(immune.combined) <- "integrated"

# Run the standard workflow for visualization and clustering
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)

# Visualization
p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2

为了并排可视化这两个条件，我们可以使用参数按群着色每个条件。

DimPlot(immune.combined, reduction = "umap", split.by = "stim")

识别保守的细胞类型标记

为了识别在各条件下保守的细胞类型标记基因，我们提供该功能。该功能为每个数据集执行基因表达检测，并使用 MetaDE R 包中的元分析方法组合 p 值。例如，我们可以计算出无论刺激条件如何，第6组（NK细胞）中保守标记的基因。

# For performing differential expression after integration, we switch back to the original data
DefaultAssay(immune.combined) <- "RNA"
nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 6, grouping.var = "stim", verbose = FALSE)
head(nk.markers)
##        CTRL_p_val CTRL_avg_log2FC CTRL_pct.1 CTRL_pct.2 CTRL_p_val_adj
## GNLY            0        6.006422      0.944      0.045              0
## FGFBP2          0        3.223246      0.503      0.020              0
## CLIC3           0        3.466418      0.599      0.024              0
## PRF1            0        2.654683      0.424      0.017              0
## CTSW            0        2.991829      0.533      0.029              0
## KLRD1           0        2.781453      0.497      0.019              0
##           STIM_p_val STIM_avg_log2FC STIM_pct.1 STIM_pct.2 STIM_p_val_adj
## GNLY    0.000000e+00        5.853573      0.956      0.060   0.000000e+00
## FGFBP2 7.275492e-161        2.200379      0.260      0.016  1.022425e-156
## CLIC3   0.000000e+00        3.549919      0.627      0.031   0.000000e+00
## PRF1    0.000000e+00        4.102686      0.862      0.057   0.000000e+00
## CTSW    0.000000e+00        3.139620      0.596      0.035   0.000000e+00
## KLRD1   0.000000e+00        2.880055      0.558      0.027   0.000000e+00
##             max_pval minimump_p_val
## GNLY    0.000000e+00              0
## FGFBP2 7.275492e-161              0
## CLIC3   0.000000e+00              0
## PRF1    0.000000e+00              0
## CTSW    0.000000e+00              0
## KLRD1   0.000000e+00              0

我们可以探索每个群的这些标记基因，并使用它们来注释我们的簇为特定的细胞类型。

FeaturePlot(immune.combined, features = c("CD3D", "SELL", "CREM", "CD8A", "GNLY", "CD79A", "FCGR3A", 
    "CCL2", "PPBP"), min.cutoff = "q9")

immune.combined <- RenameIdents(immune.combined, `0` = "CD14 Mono", `1` = "CD4 Naive T", `2` = "CD4 Memory T", 
    `3` = "CD16 Mono", `4` = "B", `5` = "CD8 T", `6` = "NK", `7` = "T activated", `8` = "DC", `9` = "B Activated", 
    `10` = "Mk", `11` = "pDC", `12` = "Eryth", `13` = "Mono/Mk Doublets", `14` = "HSPC")
DimPlot(immune.combined, label = TRUE)

参数的功能可用于查看跨条件下的保守细胞类型标记，显示表示任何给定基因在群中的表达水平和细胞百分比。在这里，我们为14个cluster中的每一个绘制2-3个强标记基因。

Idents(immune.combined) <- factor(Idents(immune.combined), levels = c("HSPC", "Mono/Mk Doublets", 
    "pDC", "Eryth", "Mk", "DC", "CD14 Mono", "CD16 Mono", "B Activated", "B", "CD8 T", "NK", "T activated", 
    "CD4 Naive T", "CD4 Memory T"))
markers.to.plot <- c("CD3D", "CREM", "HSPH1", "SELL", "GIMAP5", "CACYBP", "GNLY", "NKG7", "CCL5", 
    "CD8A", "MS4A1", "CD79A", "MIR155HG", "NME1", "FCGR3A", "VMO1", "CCL2", "S100A9", "HLA-DQA1", 
    "GPR183", "PPBP", "GNG11", "HBA2", "HBB", "TSPAN13", "IL3RA", "IGJ", "PRSS57")
DotPlot(immune.combined, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8, split.by = "stim") + 
    RotatedAxis()

识别不同条件的差异基因

现在，我们已经对齐了刺激和对照细胞，我们可以开始做比较分析，看看刺激引起的差异。观察这些变化的一种方法是绘制受刺激细胞和对照细胞的平均表达，并寻找在散点图上离群的基因。在这里，我们采取刺激和对照组的幼稚T细胞和CD14+单核细胞群的平均表达，并产生散点图，突出显示干扰素刺激剧烈反应的基因。

library(ggplot2)
library(cowplot)
theme_set(theme_cowplot())
t.cells <- subset(immune.combined, idents = "CD4 Naive T")
Idents(t.cells) <- "stim"
avg.t.cells <- as.data.frame(log1p(AverageExpression(t.cells, verbose = FALSE)$RNA))
avg.t.cells$gene <- rownames(avg.t.cells)

cd14.mono <- subset(immune.combined, idents = "CD14 Mono")
Idents(cd14.mono) <- "stim"
avg.cd14.mono <- as.data.frame(log1p(AverageExpression(cd14.mono, verbose = FALSE)$RNA))
avg.cd14.mono$gene <- rownames(avg.cd14.mono)

genes.to.label = c("ISG15", "LY6E", "IFI6", "ISG20", "MX1", "IFIT2", "IFIT1", "CXCL10", "CCL8")
p1 <- ggplot(avg.t.cells, aes(CTRL, STIM)) + geom_point() + ggtitle("CD4 Naive T Cells")
p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p2 <- ggplot(avg.cd14.mono, aes(CTRL, STIM)) + geom_point() + ggtitle("CD14 Monocytes")
p2 <- LabelPoints(plot = p2, points = genes.to.label, repel = TRUE)
p1 + p2

正如你所看到的，许多相同的基因在这两种细胞类型中都得到了调节，并且可能代表一个保守的干扰素反应途径。

我们有信心在不同条件下识别出常见的细胞类型，我们可以查看同类型细胞的基因在不同条件下有什么变化。首先，我们在meta.data 插槽中创建一个列，以同时保存细胞类型和刺激信息，并将当前标识切换到该列。然后，我们用来寻找受刺激和对照B细胞之间的不同基因。请注意，此处显示的许多top基因与我们之前绘制的核心干扰素响应基因相同。此外，我们所看到的CXCL10等基因是单核细胞和B细胞干扰素反应特有的，在此列表中也显示出非常重要的意义。

immune.combined$celltype.stim <- paste(Idents(immune.combined), immune.combined$stim, sep = "_")
immune.combined$celltype <- Idents(immune.combined)
Idents(immune.combined) <- "celltype.stim"
b.interferon.response <- FindMarkers(immune.combined, ident.1 = "B_STIM", ident.2 = "B_CTRL", verbose = FALSE)
head(b.interferon.response, n = 15)
##                 p_val avg_log2FC pct.1 pct.2     p_val_adj
## ISG15   5.398167e-155  4.5889194 0.998 0.240 7.586044e-151
## IFIT3   2.209577e-150  4.5032297 0.964 0.052 3.105118e-146
## IFI6    7.060888e-150  4.2375542 0.969 0.080 9.922666e-146
## ISG20   7.147214e-146  2.9387415 1.000 0.672 1.004398e-141
## IFIT1   7.650201e-138  4.1295888 0.914 0.032 1.075083e-133
## MX1     1.124186e-120  3.2883709 0.905 0.115 1.579819e-116
## LY6E    2.504364e-118  3.1297866 0.900 0.152 3.519383e-114
## TNFSF10 9.454398e-110  3.7783774 0.791 0.025 1.328627e-105
## IFIT2   1.672384e-105  3.6569980 0.783 0.035 2.350201e-101
## B2M      5.564362e-96  0.6100242 1.000 1.000  7.819599e-92
## PLSCR1   1.128239e-93  2.8205802 0.796 0.117  1.585514e-89
## IRF7     6.602529e-92  2.5832239 0.838 0.190  9.278534e-88
## CXCL10   4.402118e-82  5.2406913 0.639 0.010  6.186297e-78
## UBE2L6   2.995453e-81  2.1487435 0.852 0.300  4.209510e-77
## PSMB9    1.755809e-76  1.6379482 0.940 0.573  2.467438e-72

可视化基因表达变化的另一个有用的方法。显示给定基因列表的特征图，按分组变量（此处为刺激对照）进行拆分。CD3D和GNLY等基因是细胞类型标记（用于T细胞和NK/CD8 T细胞），它们几乎不受干扰素刺激的影响，并在对照和刺激组中显示类似的基因表达模式。另一方面，IFI6和ISG15是核心干扰素反应基因，在所有细胞类型中都受到相应的调节。最后，CD14和CXCL10是显示细胞类型特定干扰素反应的基因。CD14表达在CD14单核细胞刺激后减少，这可能导致监督分析框架中的误分类，强调整合分析的价值。CXCL10显示干扰素刺激后单核细胞和B细胞的明显上升调节，但其他细胞类型不显示。

FeaturePlot(immune.combined, features = c("CD3D", "GNLY", "IFI6"), split.by = "stim", max.cutoff = 3, 
    cols = c("grey", "red"))

plots <- VlnPlot(immune.combined, features = c("LYZ", "ISG15", "CXCL10"), split.by = "stim", group.by = "celltype", 
    pt.size = 0, combine = FALSE)
wrap_plots(plots = plots, ncol = 1)

使用sctransform整合数据

介绍一个改进的整合方法，该方法被命名为"sctransform"。

下面，我们演示如何修改 Seurat 整合工作流，以实现用sctransform工作流标准化数据集。命令大致相似，有几个关键差异：

• 通过SCTransform()单独（而不是在整合之前）实现数据集标准化
• 我们使用 3，000 个或更多基因来进行sctransform 下游分析。
• 在识别锚点之前运行该功能PrepSCTIntegration()
• 运行FindIntegrationAnchors()和IntegrateData()时，设置参数normalization.method为SCT
• 当运行基于 sctransform 的工作流（包括整合）时，不运行该功能ScaleData()

LoadData("ifnb")
ifnb.list <- SplitObject(ifnb, split.by = "stim")
ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform)
features <- SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 3000)
ifnb.list <- PrepSCTIntegration(object.list = ifnb.list, anchor.features = features)
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT", 
    anchor.features = features)
immune.combined.sct <- IntegrateData(anchorset = immune.anchors, normalization.method = "SCT")

immune.combined.sct <- RunPCA(immune.combined.sct, verbose = FALSE)
immune.combined.sct <- RunUMAP(immune.combined.sct, reduction = "pca", dims = 1:30)

p1 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "seurat_annotations", label = TRUE, 
    repel = TRUE)
p1 + p2

现在，数据集已经整合，您可以按照此教程中的先前步骤识别细胞类型和特定于细胞类型的响应。

文中链接

[1]

ctrl或干扰素刺激状态: https://www.nature.com/articles/nbt.4042

[2]

SeuratData: https://github.com/satijalab/seurat-data

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注